%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%         BAILOUT SIMULATIONS       %%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%            SEPTEMBER 2023         %%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear all;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Simulations varying outside assets of periphery %%
%%%                  Figure 8                     %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


Nperi = 4; 
Ncore = 6; 
N = Ncore+Ncore*Nperi;
pperiphery = 0:0.05:1;
MM = length(pperiphery); 

error = 1e-10; 

simtot = 200; 
costbiggestfirst = -999*ones(simtot, MM); 
costgreedycore = -999*ones(simtot, MM); 
costperiphery = -999*ones(simtot, MM); 
costopt = -999*ones(simtot, MM); 

for mm = 1:MM
    
    for sim = 1:simtot
        D = zeros(N,N);
        D(1:Ncore, 1:Ncore) = rand(Ncore,Ncore); % drawing intercore liabilities
        for ii = 1:Ncore % for each core bank, drawing liabilities with peripheral banks
            D(ii,Ncore+Nperi*(ii-1)+1:Ncore+Nperi*(ii-1)+Nperi) = rand(1,Nperi);
            D(Ncore+Nperi*(ii-1)+1:Ncore+Nperi*(ii-1)+Nperi,ii) = max(0,D(ii,Ncore+Nperi*(ii-1)+1:Ncore+Nperi*(ii-1)+Nperi)'-pperiphery(mm))+(1-max(0,D(ii,Ncore+Nperi*(ii-1)+1:Ncore+Nperi*(ii-1)+Nperi)'-pperiphery(mm))).*rand(Nperi,1);
        end    
        D = -(eye(size(D))-ones(size(D))).*D; %setting liabilities on self (diagonal entries) to zero
        
        DA = D*ones(N,1);
        DL = ones(1,N)*D; DL = DL'; 
        
        p = rand(N,1);
        p(Ncore+1:end) = pperiphery(mm);
        
        t0 = max(0,DL-p-DA); % transfers to guarantee weak balance which are made first
        
        % computing worst equilibrium absent additional transfers
        S0 = zeros(N,1); % set of solvent banks
        S0bis = p+t0>=DL; 
        while sum(S0bis-S0)>0
            S0 = S0bis; 
            S0bis = p+t0+D*S0>=DL;
        end   

        % biggest banks first
        sizeasset = p+DA; 
        tbiggestfirst = zeros(N,1); 

        S=S0; 
        while sum(S)<N
            BC = DL-p -t0-tbiggestfirst-D*S;
            ratio = zeros(N,1);
            ratio(S==0) = sizeasset(S==0);
            [~,idx] = max(ratio);
            tbiggestfirst(idx) = BC(idx);
            Sbis = p+t0+tbiggestfirst+D*S>=DL-error;
            while max(Sbis-S)>0
                S=Sbis; 
                Sbis = p+t0+tbiggestfirst+D*S>=DL-error;
            end    
        end
        costbiggestfirst(sim,mm) = sum(tbiggestfirst); 
        
        % greedy on core
        IBV = -999*ones(N,1); 
        tgreedycore = zeros(N,1); 
        orderbailout = zeros(Ncore,1);

        S=S0; 
        while sum(S)<N
            for ii = 1:Ncore
                IBV(ii) = sum(min(D(:,ii),max(0,DL-p -t0-tgreedycore-(D*S)))); 
            end
            BC = DL-p -t0-tgreedycore-D*S;
            ratio = zeros(N,1);
            ratio(S==0) = IBV(S==0)./BC(S==0);
            ratio(Ncore+1:end)=0;
            [~,idx] = max(ratio);
            tgreedycore(idx) = BC(idx);
            orderbailout(sum(orderbailout>0)+1)=idx;
            Sbis = p+t0+tgreedycore+D*S>=DL-error;
            while max(Sbis-S)>0
                S=Sbis; 
                Sbis = p+t0+tgreedycore+D*S>=DL-error;
            end    
        end
        costgreedycore(sim,mm) = sum(tgreedycore); 
        
        % opt periphery before each core bank

        IBV = -999*ones(N,1); 
        tperiphery = zeros(N,1); 
        count = 0; 

        S=S0; 
        while sum(S)<N
            BC = max(DL-p -t0-tperiphery-D*S,0);
            count = count+1;
            idx = orderbailout(count);
            tperiphery(idx) = BC(idx);
            for nn = 1:Nperi
                tempbailout = nchoosek(1:Nperi,nn);
                tempcost = -999*ones(size(tempbailout,1),1);
                for kk = 1:size(tempbailout,1)
                    tempcost(kk) = sum(BC(Ncore + (idx-1)*Nperi+tempbailout(kk,:)'));
                    tempcost(kk) = tempcost(kk) + max(0,BC(idx) - sum(D(idx,Ncore + (idx-1)*Nperi+tempbailout(kk,:)).*(ones(1,nn)-S(Ncore + (idx-1)*Nperi+tempbailout(kk,:))')));
                end    
                if min(tempcost)<tperiphery(idx) + sum(tperiphery(Ncore + (idx-1)*Nperi +1 :Ncore + (idx-1)*Nperi +Nperi))-error
                    [~,lowidx] = min(tempcost);
                    tperiphery(Ncore + (idx-1)*Nperi+tempbailout(lowidx,:)') = BC(Ncore + (idx-1)*Nperi+tempbailout(lowidx,:)'); 
                    tperiphery(idx) =max(0,BC(idx) - sum(D(idx,Ncore + (idx-1)*Nperi+tempbailout(lowidx,:)).*(ones(1,nn)-S(Ncore + (idx-1)*Nperi+tempbailout(lowidx,:))')));
                end 
            end
            Sbis = p+t0+tperiphery+D*S>=DL-error;
            while max(Sbis-S)>0
                S=Sbis; 
                Sbis = p+t0+tperiphery+D*S>=DL-error;
            end  
        end
        costperiphery(sim,mm) = sum(tperiphery); 

        % OPT solution

        topt = 999*ones(N,1); 

        policies = perms(1:Ncore);
        PP = size(policies,1); 
        for pp=1:PP
          count = 0;
          S=S0; 
          toptbis = zeros(N,1); 
          while sum(S)<N
                BC = max(DL-p -t0-toptbis-D*S,0);
                count = count+1;
                idx = policies(pp,count);
                toptbis(idx) = BC(idx);
                for nn = 1:Nperi
                    tempbailout = nchoosek(1:Nperi,nn);
                    tempcost = -999*ones(size(tempbailout,1),1);
                    for kk = 1:size(tempbailout,1)
                        tempcost(kk) = sum(BC(Ncore + (idx-1)*Nperi+tempbailout(kk,:)'));
                        tempcost(kk) = tempcost(kk) + max(0,BC(idx) - sum(D(idx,Ncore + (idx-1)*Nperi+tempbailout(kk,:)).*(ones(1,nn)-S(Ncore + (idx-1)*Nperi+tempbailout(kk,:))')));
                    end    
                    if min(tempcost)<toptbis(idx) + sum(toptbis(Ncore + (idx-1)*Nperi +1 :Ncore + (idx-1)*Nperi +Nperi))-error
                        [~,lowidx] = min(tempcost);
                        toptbis(Ncore + (idx-1)*Nperi+tempbailout(lowidx,:)') = BC(Ncore + (idx-1)*Nperi+tempbailout(lowidx,:)'); 
                        toptbis(idx) =max(0,BC(idx) - sum(D(idx,Ncore + (idx-1)*Nperi+tempbailout(lowidx,:)).*(ones(1,nn)-S(Ncore + (idx-1)*Nperi+tempbailout(lowidx,:))')));
                    end 
                end
                Sbis = p+t0+toptbis+D*S>=DL-error;
                while max(Sbis-S)>0
                    S=Sbis; 
                    Sbis = p+t0+toptbis+D*S>=DL-error;
                end  
            end
            if sum(toptbis)<sum(topt)-error
                topt = toptbis;
            end    
        end    
        costopt(sim,mm) = sum(topt); 

    end
end    

%% Plots %%

figure('DefaultAxesFontSize',10); 
plot(pperiphery, mean(costbiggestfirst),'--square', 'Color',1/255*[255, 173, 168],'LineWidth',2.5); hold on; 
plot(pperiphery, mean(costgreedycore),'--v', 'Color',1/255*[235, 119, 116],'LineWidth',2.5); hold on; 
plot(pperiphery, mean(costperiphery),'--^','Color',1/255*[146, 39, 46],'LineWidth',2.5); hold on; 
plot(pperiphery, mean(costopt),'--o','Color',1/255*[61, 0, 0],'LineWidth',2.5); hold on; 
legend('Core First, by size', 'Core First, by $\frac{IBV}{BC}$','Periphery First', 'Optimum', 'Interpreter','latex', 'location', 'northeast','FontSize', 17);
xlabel('Assets of Peripheral Banks', 'Interpreter','latex','FontSize', 19);
ylabel('Total Bailout Costs', 'Interpreter','latex','FontSize', 19);
title('Recovery Rate $(1-a)=0$', 'Interpreter','latex','FontSize', 19);


% figure('DefaultAxesFontSize',12); 
% plot(pperiphery, mean(costgreedycore),'--square', 'Color',1/255*[255, 145, 142],'LineWidth',2.5); hold on; 
% plot(pperiphery, mean(costperiphery),'--^','Color',1/255*[175, 66, 68],'LineWidth',2.5); hold on; 
% plot(pperiphery, mean(costopt),'--o','Color',1/255*[61, 0, 0],'LineWidth',2.5); hold on; 
% legend('Core First','Periphery First', 'Optimum', 'Interpreter','latex', 'location', 'northeast','FontSize', 20);
% xlabel('Assets of Peripheral Banks', 'Interpreter','latex','FontSize', 20);
% ylabel('Total Bailout Costs', 'Interpreter','latex','FontSize', 20);
% title('$n_C=6$, $n_P=4$','Interpreter','latex','FontSize', 20);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%   Same as above but with partial recovery rates  %%
%%% Simulations varying outside assets of periphery %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


Nperi = 4; 
Ncore = 6; 
N = Ncore+Ncore*Nperi;
costpara_a = 0.25; %recovery rate
pperiphery = [0:0.05:0.40 0.55:0.15:1];
MM = length(pperiphery); 

error = 1e-10; 

simtot = 200; 
costbiggestfirst  = -999*ones(simtot, MM); 
costgreedycore = -999*ones(simtot, MM); 
costperiphery = -999*ones(simtot, MM); 
costopt = -999*ones(simtot, MM); 

for mm = 1:MM

    for sim = 201:500
        D = zeros(N,N);
        D(1:Ncore, 1:Ncore) = rand(Ncore,Ncore); % drawing intercore liabilities
        for ii = 1:Ncore % for each core bank, drawing liabilities with peripheral banks
            D(ii,Ncore+Nperi*(ii-1)+1:Ncore+Nperi*(ii-1)+Nperi) = rand(1,Nperi);
            D(Ncore+Nperi*(ii-1)+1:Ncore+Nperi*(ii-1)+Nperi,ii) = max(0,D(ii,Ncore+Nperi*(ii-1)+1:Ncore+Nperi*(ii-1)+Nperi)'-pperiphery(mm))+(1-max(0,D(ii,Ncore+Nperi*(ii-1)+1:Ncore+Nperi*(ii-1)+Nperi)'-pperiphery(mm))).*rand(Nperi,1);
        end    
        D = -(eye(size(D))-ones(size(D))).*D; %setting liabilities on self (diagonal entries) to zero
        
        DA = D*ones(N,1);
        DL = ones(1,N)*D; DL = DL'; 
        Pi = D./(repmat(DL',N,1)); 

        p = rand(N,1);
        p(Ncore+1:end) = pperiphery(mm);
        
        t0 = max(0,DL-p-DA); % transfers to guarantee weak balance which are made first
        
        % computing worst equilibrium absent additional transfers
        paymentsvec = zeros(1,N);
        payments = Pi.*repmat(paymentsvec,N,1);
        S0 = p+t0+sum(payments,2)>=DL-error; 

        paymentsvecbis = zeros(1,N);
        paymentsvecbis(S0==1) = DL(S0==1);
        tempassets = p+t0+sum(payments,2);
        paymentsvecbis(S0==0) = (1-costpara_a)*tempassets(S0==0);
        while max(abs(paymentsvecbis-paymentsvec))>error
            paymentsvec = paymentsvecbis; 
            payments = Pi.*repmat(paymentsvec,N,1);
            S0 = p+t0+sum(payments,2)>=DL-error; 
            paymentsvecbis(S0==1) = DL(S0==1);
            tempassets = p+t0+sum(payments,2);
            paymentsvecbis(S0==0) = (1-costpara_a)*tempassets(S0==0);
        end          
        paymentsworsteq = payments;

        % biggest banks first
        sizeasset = p+DA; 
        tbiggestfirst = zeros(N,1); 

        S=S0; 
        while sum(S)<N
            BC = DL-p -t0-tbiggestfirst-sum(payments,2);
            ratio = zeros(N,1);
            ratio(S==0) = sizeasset(S==0);
            [~,idx] = max(ratio);
            tbiggestfirst(idx) = BC(idx);

            paymentsvec = zeros(1,N);
            payments = Pi.*repmat(paymentsvec,N,1);
            S = p+t0+tbiggestfirst+sum(payments,2)>=DL-error; 
    
            paymentsvecbis = zeros(1,N);
            paymentsvecbis(S==1) = DL(S==1);
            tempassets = p+t0+tbiggestfirst+sum(payments,2);
            paymentsvecbis(S==0) = (1-costpara_a)*tempassets(S==0);
            while max(abs(paymentsvecbis-paymentsvec))>error
                paymentsvec = paymentsvecbis; 
                payments = Pi.*repmat(paymentsvec,N,1);
                S = p+t0+tbiggestfirst+sum(payments,2)>=DL-error; 
                paymentsvecbis(S==1) = DL(S==1);
                tempassets = p+t0+tbiggestfirst+sum(payments,2);
                paymentsvecbis(S==0) = (1-costpara_a)*tempassets(S==0);
            end           
        end
        costbiggestfirst(sim,mm) = sum(tbiggestfirst); 
        
        % greedy on core
        IBV = -999*ones(N,1); 
        tgreedycore = zeros(N,1); 
        orderbailout = zeros(Ncore,1);

        S=S0; 
        payments = paymentsworsteq;
        while sum(S)<N
            for ii = 1:Ncore
                IBV(ii) = sum(min(D(:,ii),max(0,DL-p -t0-tgreedycore-sum(payments,2)))); 
            end
            BC = DL-p -t0-tgreedycore-sum(payments,2);
            ratio = zeros(N,1);
            ratio(S==0) = IBV(S==0)./BC(S==0);
            ratio(Ncore+1:end)=0;
            [~,idx] = max(ratio);
            tgreedycore(idx) = BC(idx);
            orderbailout(sum(orderbailout>0)+1)=idx;

            paymentsvec = zeros(1,N);
            payments = Pi.*repmat(paymentsvec,N,1);
            S = p+t0+tgreedycore+sum(payments,2)>=DL-error; 
    
            paymentsvecbis = zeros(1,N);
            paymentsvecbis(S==1) = DL(S==1);
            tempassets = p+t0+tgreedycore+sum(payments,2);
            paymentsvecbis(S==0) = (1-costpara_a)*tempassets(S==0);
            while max(abs(paymentsvecbis-paymentsvec))>error
                paymentsvec = paymentsvecbis; 
                payments = Pi.*repmat(paymentsvec,N,1);
                S = p+t0+tgreedycore+sum(payments,2)>=DL-error; 
                paymentsvecbis(S==1) = DL(S==1);
                tempassets = p+t0+tgreedycore+sum(payments,2);
                paymentsvecbis(S==0) = (1-costpara_a)*tempassets(S==0);
            end           
        end
        costgreedycore(sim,mm) = sum(tgreedycore); 
        
        % opt periphery before each core bank

        IBV = -999*ones(N,1); 
        tperiphery = zeros(N,1); 
        count = 0; 

        S=S0; 
        payments = paymentsworsteq;

        while sum(S)<N
            BC = max(DL-p -t0-tperiphery-sum(payments,2),0);
            count = count+1;
            idx = orderbailout(count);
            tperiphery(idx) = BC(idx);
            for nn = 1:Nperi
                tempbailout = nchoosek(1:Nperi,nn);
                tempcost = -999*ones(size(tempbailout,1),1);
                for kk = 1:size(tempbailout,1)
                    tempcost(kk) = sum(BC(Ncore + (idx-1)*Nperi+tempbailout(kk,:)'));
                    tempcost(kk) = tempcost(kk) + max(0,BC(idx) - sum(D(idx,Ncore + (idx-1)*Nperi+tempbailout(kk,:))-payments(idx,Ncore + (idx-1)*Nperi+tempbailout(kk,:))));
                end    
                if min(tempcost)<tperiphery(idx) + sum(tperiphery(Ncore + (idx-1)*Nperi +1 :Ncore + (idx-1)*Nperi +Nperi))-error
                    [~,lowidx] = min(tempcost);
                    tperiphery(Ncore + (idx-1)*Nperi+tempbailout(lowidx,:)') = BC(Ncore + (idx-1)*Nperi+tempbailout(lowidx,:)'); 
                    tperiphery(idx) =max(0,BC(idx) - sum(D(idx,Ncore + (idx-1)*Nperi+tempbailout(lowidx,:))-payments(idx,Ncore + (idx-1)*Nperi+tempbailout(lowidx,:))));
                end 
            end
            paymentsvec = zeros(1,N);
            payments = Pi.*repmat(paymentsvec,N,1);
            S = p+t0+tperiphery+sum(payments,2)>=DL-error; 
    
            paymentsvecbis = zeros(1,N);
            paymentsvecbis(S==1) = DL(S==1);
            tempassets = p+t0+tperiphery+sum(payments,2);
            paymentsvecbis(S==0) = (1-costpara_a)*tempassets(S==0);
            while max(abs(paymentsvecbis-paymentsvec))>error
                paymentsvec = paymentsvecbis; 
                payments = Pi.*repmat(paymentsvec,N,1);
                S = p+t0+tperiphery+sum(payments,2)>=DL-error; 
                paymentsvecbis(S==1) = DL(S==1);
                tempassets = p+t0+tperiphery+sum(payments,2);
                paymentsvecbis(S==0) = (1-costpara_a)*tempassets(S==0);
            end  
        end
        costperiphery(sim,mm) = sum(tperiphery); 

        % OPT solution

        topt = 999*ones(N,1); 

        policies = perms(1:Ncore);
        PP = size(policies,1); 
        for pp=1:PP
          count = 0;
          S=S0; 
          payments = paymentsworsteq;
          toptbis = zeros(N,1); 
          while sum(S)<N
                BC = max(DL-p -t0-toptbis-sum(payments,2),0);
                count = count+1;
                idx = policies(pp,count);
                toptbis(idx) = BC(idx);
                for nn = 1:Nperi
                    tempbailout = nchoosek(1:Nperi,nn);
                    tempcost = -999*ones(size(tempbailout,1),1);
                    for kk = 1:size(tempbailout,1)
                        tempcost(kk) = sum(BC(Ncore + (idx-1)*Nperi+tempbailout(kk,:)'));
                        tempcost(kk) = tempcost(kk) + max(0,BC(idx) - sum(D(idx,Ncore + (idx-1)*Nperi+tempbailout(kk,:))-payments(idx,Ncore + (idx-1)*Nperi+tempbailout(kk,:))));
                    end    
                    if min(tempcost)<toptbis(idx) + sum(toptbis(Ncore + (idx-1)*Nperi +1 :Ncore + (idx-1)*Nperi +Nperi))-error
                        [~,lowidx] = min(tempcost);
                        toptbis(Ncore + (idx-1)*Nperi+tempbailout(lowidx,:)') = BC(Ncore + (idx-1)*Nperi+tempbailout(lowidx,:)'); 
                        toptbis(idx) =max(0,BC(idx) - sum(D(idx,Ncore + (idx-1)*Nperi+tempbailout(lowidx,:))-payments(idx,Ncore + (idx-1)*Nperi+tempbailout(lowidx,:))));
                    end 
                end
 
                paymentsvec = zeros(1,N);
                payments = Pi.*repmat(paymentsvec,N,1);
                S = p+t0+toptbis+sum(payments,2)>=DL-error; 
        
                paymentsvecbis = zeros(1,N);
                paymentsvecbis(S==1) = DL(S==1);
                tempassets = p+t0+toptbis+sum(payments,2);
                paymentsvecbis(S==0) = (1-costpara_a)*tempassets(S==0);
                while max(abs(paymentsvecbis-paymentsvec))>error
                    paymentsvec = paymentsvecbis; 
                    payments = Pi.*repmat(paymentsvec,N,1);
                    S = p+t0+toptbis+sum(payments,2)>=DL-error; 
                    paymentsvecbis(S==1) = DL(S==1);
                    tempassets = p+t0+toptbis+sum(payments,2);
                    paymentsvecbis(S==0) = (1-costpara_a)*tempassets(S==0);
                end  
        end
            if sum(toptbis)<sum(topt)-error
                topt = toptbis;
            end    
        end    
        costopt(sim,mm) = sum(topt); 
    end
end    

%% Plots %%

figure('DefaultAxesFontSize',10); 
plot(pperiphery, mean(costbiggestfirst(:,1:end)),'--square', 'Color',1/255*[255, 173, 168],'LineWidth',2.5); hold on; 
plot(pperiphery, mean(costgreedycore(:,1:end)),'--v', 'Color',1/255*[235, 119, 116],'LineWidth',2.5); hold on; 
plot(pperiphery, mean(costperiphery(:,1:end)),'--^','Color',1/255*[146, 39, 46],'LineWidth',2.5); hold on; 
plot(pperiphery, mean(costopt(:,1:end)),'--o','Color',1/255*[61, 0, 0],'LineWidth',2.5); hold on; 
legend('Core First, by Size', 'Core First, by $\frac{IBV}{BC}$','Periphery First', 'Optimum', 'Interpreter','latex', 'location', 'northeast','FontSize', 17);
xlabel('Assets of Peripheral Banks', 'Interpreter','latex','FontSize', 19);
ylabel('Total Bailout Costs', 'Interpreter','latex','FontSize', 19);
title('Recovery Rate $(1-a)=0.75$', 'Interpreter','latex','FontSize', 19);
%ylim([0,9]);%xlim([1,MM]);   yticks(0.5:0.1:1); 



%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%% Simulations Varying Asymmetry of Core Banks %%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%                  Figure 9                   %%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


Nperi = 4; 
Ncore = 6; 
N = Ncore+Ncore*Nperi;
costpara_a = 0.25; 
asymmetry_vec = 0:0.1:1; 
MM = length(asymmetry_vec); 
pperiphery = 0.1; 

error = 1e-10; 

simtot = 1000; 
costperiphery = -999*ones(simtot, MM); 
costopt = -999*ones(simtot, MM); 

for mm = 1:MM
    % drawing the network 
    for sim = 1:simtot
        D = zeros(N,N);
        a = asymmetry_vec(mm);
        temp = rand(1,1)*ones(Ncore,Ncore); 
        D(1:Ncore, 1:Ncore) = (1-a)*temp+a*rand(Ncore,Ncore); % drawing intercore liabilities
        D(1,Ncore+Nperi*(1-1)+1:Ncore+Nperi*(1-1)+Nperi) = rand(1,Nperi);
        D(Ncore+Nperi*(1-1)+1:Ncore+Nperi*(1-1)+Nperi,1) = rand(Nperi,1);
        D(Ncore+Nperi*(1-1)+1:Ncore+Nperi*(1-1)+Nperi,1) = max(0,D(1,Ncore+Nperi*(1-1)+1:Ncore+Nperi*(1-1)+Nperi)'-pperiphery)+(1-max(0,D(1,Ncore+Nperi*(1-1)+1:Ncore+Nperi*(1-1)+Nperi)'-pperiphery)).*rand(Nperi,1);
        for ii = 2:Ncore % for each core bank, drawing liabilities with peripheral banks
            D(ii,Ncore+Nperi*(ii-1)+1:Ncore+Nperi*(ii-1)+Nperi) = (1-a)*D(1,Ncore+Nperi*(1-1)+1:Ncore+Nperi*(1-1)+Nperi)+a*rand(1,Nperi);
            D(Ncore+Nperi*(ii-1)+1:Ncore+Nperi*(ii-1)+Nperi,ii) = (1-a)*D(Ncore+Nperi*(1-1)+1:Ncore+Nperi*(1-1)+Nperi,1)+a*rand(Nperi,1);
            D(Ncore+Nperi*(ii-1)+1:Ncore+Nperi*(ii-1)+Nperi,ii) = (1-a)*D(Ncore+Nperi*(1-1)+1:Ncore+Nperi*(1-1)+Nperi,1)+a*max(0,D(ii,Ncore+Nperi*(ii-1)+1:Ncore+Nperi*(ii-1)+Nperi)'-pperiphery)+(1-max(0,D(ii,Ncore+Nperi*(ii-1)+1:Ncore+Nperi*(ii-1)+Nperi)'-pperiphery)).*rand(Nperi,1);
        end    
        D = -(eye(size(D))-ones(size(D))).*D; %setting liabilities on self (diagonal entries) to zero
        
        DA = D*ones(N,1);
        DL = ones(1,N)*D; DL = DL'; 
        Pi = D./(repmat(DL',N,1)); 
        
        p = (1-a)*rand(1,1)+a*rand(N,1);
        pubarperi = 0.2;
        plbarperi = 0; 
        temp = repmat(unifrnd(plbarperi,pubarperi,Nperi,1),Ncore,1);
        p(Ncore+1:end) = (1-a)*temp+a*unifrnd(plbarperi,pubarperi,Ncore*Nperi,1);
        p(Ncore+1:end) = pperiphery;

        t0 = max(0,DL-p-DA); % transfers to guarantee weak balance which are made first
        
        % computing worst equilibrium absent additional transfers
        paymentsvec = zeros(1,N);
        payments = Pi.*repmat(paymentsvec,N,1);
        S0 = p+t0+sum(payments,2)>=DL-error; 

        paymentsvecbis = zeros(1,N);
        paymentsvecbis(S0==1) = DL(S0==1);
        tempassets = p+t0+sum(payments,2);
        paymentsvecbis(S0==0) = (1-costpara_a)*tempassets(S0==0);
        while max(abs(paymentsvecbis-paymentsvec))>error
            paymentsvec = paymentsvecbis; 
            payments = Pi.*repmat(paymentsvec,N,1);
            S0 = p+t0+sum(payments,2)>=DL-error; 
            paymentsvecbis(S0==1) = DL(S0==1);
            tempassets = p+t0+sum(payments,2);
            paymentsvecbis(S0==0) = (1-costpara_a)*tempassets(S0==0);
        end          
        paymentsworsteq = payments;

        % greedy on core
        IBV = -999*ones(N,1); 
        tgreedycore = zeros(N,1); 
        orderbailout = zeros(Ncore,1);

        S=S0; 
        payments = paymentsworsteq;
        while sum(S)<N
            for ii = 1:Ncore
                IBV(ii) = sum(min(D(:,ii),max(0,DL-p -t0-tgreedycore-sum(payments,2)))); 
            end
            BC = DL-p -t0-tgreedycore-sum(payments,2);
            ratio = zeros(N,1);
            ratio(S==0) = IBV(S==0)./BC(S==0);
            ratio(Ncore+1:end)=0;
            [~,idx] = max(ratio);
            tgreedycore(idx) = BC(idx);
            orderbailout(sum(orderbailout>0)+1)=idx;

            paymentsvec = zeros(1,N);
            payments = Pi.*repmat(paymentsvec,N,1);
            S = p+t0+tgreedycore+sum(payments,2)>=DL-error; 
    
            paymentsvecbis = zeros(1,N);
            paymentsvecbis(S==1) = DL(S==1);
            tempassets = p+t0+tgreedycore+sum(payments,2);
            paymentsvecbis(S==0) = (1-costpara_a)*tempassets(S==0);
            while max(abs(paymentsvecbis-paymentsvec))>error
                paymentsvec = paymentsvecbis; 
                payments = Pi.*repmat(paymentsvec,N,1);
                S = p+t0+tgreedycore+sum(payments,2)>=DL-error; 
                paymentsvecbis(S==1) = DL(S==1);
                tempassets = p+t0+tgreedycore+sum(payments,2);
                paymentsvecbis(S==0) = (1-costpara_a)*tempassets(S==0);
            end           
        end
        
        % opt periphery before each core bank

        IBV = -999*ones(N,1); 
        tperiphery = zeros(N,1); 
        count = 0; 

        S=S0; 
        payments = paymentsworsteq;

        while sum(S)<N
            BC = max(DL-p -t0-tperiphery-sum(payments,2),0);
            count = count+1;
            idx = orderbailout(count);
            tperiphery(idx) = BC(idx);
            for nn = 1:Nperi
                tempbailout = nchoosek(1:Nperi,nn);
                tempcost = -999*ones(size(tempbailout,1),1);
                for kk = 1:size(tempbailout,1)
                    tempcost(kk) = sum(BC(Ncore + (idx-1)*Nperi+tempbailout(kk,:)'));
                    tempcost(kk) = tempcost(kk) + max(0,BC(idx) - sum(D(idx,Ncore + (idx-1)*Nperi+tempbailout(kk,:))-payments(idx,Ncore + (idx-1)*Nperi+tempbailout(kk,:))));
                end    
                if min(tempcost)<tperiphery(idx) + sum(tperiphery(Ncore + (idx-1)*Nperi +1 :Ncore + (idx-1)*Nperi +Nperi))-error
                    [~,lowidx] = min(tempcost);
                    tperiphery(Ncore + (idx-1)*Nperi+tempbailout(lowidx,:)') = BC(Ncore + (idx-1)*Nperi+tempbailout(lowidx,:)'); 
                    tperiphery(idx) =max(0,BC(idx) - sum(D(idx,Ncore + (idx-1)*Nperi+tempbailout(lowidx,:))-payments(idx,Ncore + (idx-1)*Nperi+tempbailout(lowidx,:))));
                end 
            end
 
            paymentsvec = zeros(1,N);
            payments = Pi.*repmat(paymentsvec,N,1);
            S = p+t0+tperiphery+sum(payments,2)>=DL-error; 
    
            paymentsvecbis = zeros(1,N);
            paymentsvecbis(S==1) = DL(S==1);
            tempassets = p+t0+tperiphery+sum(payments,2);
            paymentsvecbis(S==0) = (1-costpara_a)*tempassets(S==0);
            while max(abs(paymentsvecbis-paymentsvec))>error
                paymentsvec = paymentsvecbis; 
                payments = Pi.*repmat(paymentsvec,N,1);
                S = p+t0+tperiphery+sum(payments,2)>=DL-error; 
                paymentsvecbis(S==1) = DL(S==1);
                tempassets = p+t0+tperiphery+sum(payments,2);
                paymentsvecbis(S==0) = (1-costpara_a)*tempassets(S==0);
            end  
        end
        costperiphery(sim,mm) = sum(tperiphery); 

        % OPT solution

        topt = 999*ones(N,1); 

        policies = perms(1:Ncore);
        PP = size(policies,1); 
        for pp=1:PP
          count = 0;
          S=S0; 
          payments = paymentsworsteq;
          toptbis = zeros(N,1); 
          while sum(S)<N
                BC = max(DL-p -t0-toptbis-sum(payments,2),0);
                count = count+1;
                idx = policies(pp,count);
                toptbis(idx) = BC(idx);
                for nn = 1:Nperi
                    tempbailout = nchoosek(1:Nperi,nn);
                    tempcost = -999*ones(size(tempbailout,1),1);
                    for kk = 1:size(tempbailout,1)
                        tempcost(kk) = sum(BC(Ncore + (idx-1)*Nperi+tempbailout(kk,:)'));
                        %tempcost(kk) = tempcost(kk) + max(0,BC(idx) - sum(D(idx,Ncore + (idx-1)*Nperi+tempbailout(kk,:)).*(ones(1,nn)-S(Ncore + (idx-1)*Nperi+tempbailout(kk,:))')));
                        tempcost(kk) = tempcost(kk) + max(0,BC(idx) - sum(D(idx,Ncore + (idx-1)*Nperi+tempbailout(kk,:))-payments(idx,Ncore + (idx-1)*Nperi+tempbailout(kk,:))));
                    end    
                    if min(tempcost)<toptbis(idx) + sum(toptbis(Ncore + (idx-1)*Nperi +1 :Ncore + (idx-1)*Nperi +Nperi))-error
                        [~,lowidx] = min(tempcost);
                        toptbis(Ncore + (idx-1)*Nperi+tempbailout(lowidx,:)') = BC(Ncore + (idx-1)*Nperi+tempbailout(lowidx,:)'); 
                        toptbis(idx) =max(0,BC(idx) - sum(D(idx,Ncore + (idx-1)*Nperi+tempbailout(lowidx,:))-payments(idx,Ncore + (idx-1)*Nperi+tempbailout(lowidx,:))));
                    end 
                end

                paymentsvec = zeros(1,N);
                payments = Pi.*repmat(paymentsvec,N,1);
                S = p+t0+toptbis+sum(payments,2)>=DL-error; 
        
                paymentsvecbis = zeros(1,N);
                paymentsvecbis(S==1) = DL(S==1);
                tempassets = p+t0+toptbis+sum(payments,2);
                paymentsvecbis(S==0) = (1-costpara_a)*tempassets(S==0);
                while max(abs(paymentsvecbis-paymentsvec))>error
                    paymentsvec = paymentsvecbis; 
                    payments = Pi.*repmat(paymentsvec,N,1);
                    S = p+t0+toptbis+sum(payments,2)>=DL-error; 
                    paymentsvecbis(S==1) = DL(S==1);
                    tempassets = p+t0+toptbis+sum(payments,2);
                    paymentsvecbis(S==0) = (1-costpara_a)*tempassets(S==0);
                end  
        end
            if sum(toptbis)<sum(topt)-error
                topt = toptbis;
            end    
        end    
        costopt(sim,mm) = sum(topt);
    end
end  

%% Plots %%

figure('DefaultAxesFontSize',10); 
plot(asymmetry_vec(1:end), mean(costperiphery(:,1:end)),'--^','Color',1/255*[146, 39, 46],'LineWidth',2.5); hold on; 
plot(asymmetry_vec(1:end), mean(costopt(:,1:end)),'--o','Color',1/255*[61, 0, 0],'LineWidth',2.5); hold on; 
legend('Periphery First', 'Optimum', 'Interpreter','latex', 'location', 'northeast','FontSize', 19);
xlabel('Asymmetry of Core Banks', 'Interpreter','latex','FontSize', 19);
ylabel('Total Bailout Costs', 'Interpreter','latex','FontSize', 19);
title('Recovery Rate $(1-a)=0.75$', 'Interpreter','latex','FontSize', 19);


